Spin in Density-Functional Theory 

Christoph R. JacotQ and Markus Reihei]^ 

^Karlsruhe Institute of Technology (KIT), 
Center for Functional Nanostructures and Institute of Physical Chemistry, 
Wolfgang-Gaede-StraBe la, 76131 Karlsruhe, Germany 

^ETH Zurich, Laboratorium fiir Physikalische Chemie, 
Wolfgang-Pauli-Strasse 10, 8093 Zurich, Switzerland 

Abstract 

The accurate description of open-shell molecules, in particular of transition metal 
complexes and clusters, is still an important challenge for quantum chemistry. While 
density-functional theory (DFT) is widely applied in this area, the sometimes severe 
limitations of its currently available approximate realizations often preclude its ap- 
plication as a predictive theory. Here, we review the foundations of DFT applied to 
open-shell systems, both within the nonrelativistic and the relativistic framework. 
In particular, we provide an in-depth discussion of the exact theory, with a focus 
on the role of the spin density and possibilities for targeting specific spin states. 
It turns out that different options exist for setting up Kohn-Sham DFT schemes 
for open-shell systems, which imply different definitions of the exchange-correlation 
energy functional and lead to different exact conditions on this functional. Finally, 
we suggest possible directions for future developments. 



Date: October 25, 2012 

Status: published in Int. J. Quantum Chem. 112, 3661-3684 (2012). 



DOI: 



http : //dx . doi . org/10 . 1002/qua . 24309 



^E-Mail: christopli.jacob@kit.edu 
^E-Mail: markus. reilier@phys. chem. ethz.ch 



1 



Table of Contents Graphics and Text 



SpinDFT 

p(r) 

Qir) {S^) 

In this Tutorial Review, we outline the foundations of density-functional theory (DFT) 
applied to open-shell systems, both in the non-relativistic case and within the relativistic 
theory. The role of the spin density as well as possibilities for targeting specific spin states 
are discussed, and we suggest some possible future directions for Spin-DFT. 
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1 Introduction 



Open-shell molecules such as, for example, radicals or transition metal complexes and 
clusters, feature a measurable magnetic moment that originates from their electronic 
structure. In fact, the electronic spin gives rise to a magnetic moment that makes such 
molecular systems functional for various purposes. For instance, organic radicals can be 
employed as spin probes in biomolecules [Tj and are of interest as building blocks for 
molecular spintronics devices |2|[3], single- molecule magnets have the potential to act as 
molecular qubits for quantum information processing |4| , and open-shell transition metal 
compounds serve as catalytic centers in (bio-)inorganic chemistry [5]-[8], where a change 
in the spin state can be an essential step in the catalytic cycle lol. 



Consequently, a first-principles theory that is useful for descriptive and analytic pur- 
poses and that has the potential to be a predictive tool in theoretical studies on such 
chemical systems must consider the spin properties of the electronic structure. While for 
closed-shell systems, quantum chemical methods — both wavefunction theory for accurate 



calculations on small molecules 10 and density-functional theory (DFT) for studies on 
complex chemical systems flT\ — offer such predictive tools, the situation is less satisfac- 
tory for open-shell systems, in particular for transition metal complexes and clusters fl2]. 



With wavefunction based methods, a multi-reference treatment is in general mandatory for 
open-shell systems. In particular, the complete active space self-consistent field (CASSCF) 
method, usually in combination with second-order perturbation theory (CASPT2), has 
been employed to study transition metal complexes (for examples, see Refs. jT3}jT8] ). 
However, the factorial scaling with the size of the active space puts rather severe limits 
on the size of the active space, which prevents most applications to polynuclear transition 
metal complexes and clusters. Novel approaches, such as the density matrix renormaliza- 



tion group (DMRG) algorithm 19 20 and its generalizations 21 might make it possible 



to overcome this limitation, although the molecular sizes that can be studied are clearly 
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much smaller compared to those accessible to DFT methods. 



Therefore, DFT is usually the method of choice in theoretical studies of transition-metal 
catalysis as well as molecular and spectroscopic properties of open-shell molecular systems 



22 25 . Despite much success, it has also become clear that for open-shell systems, DFT 



with the currently available approximate functionals shows a number of shortcomings. In 
addition to inaccuracies in predicting energies, geometries, and molecular properties (for 



a case study, see, e.g., Ref. 26 and for overviews, see, e.g., Refs. 11,27,28]), a severe 
limitation are unsystematic errors in the prediction of the relative energies of different spin 



states 29 -36 . Moreover, the spin density — which serves as an additional fundamental 



quantity in the spin-DFT formalism commonly employed for open-shell systems — is 



qualitatively incorrect in some cases 15,37-39 . To make things even worse, the treatment 



of low-spin states usually requires the use of a broken-symmetry description 40-43 



which provides an unphysical spin density by construction (see, e.g., Refs. 24,44 for 



a discussion). This precludes the simple prediction of spectroscopic properties depending 



on the spin density (for schemes to address this difficulty, see, e.g., Refs. 45-48]). 



Consequently, the development of better approximate DFT methods for open-shell sys- 
tems is currently still one of the most important and challenging topics in theoretical 



chemistry 12,28,49 . To make progress is this area, it is important to understand the 



exact theory underlying DFT for open-shell systems. While for the closed-shell case, ex- 



haustive presentations of this theory exists in several textbooks 50 -53 , this is not the 



case for open-shell systems, which are often only mentioned in passing in these accounts. 
Here, we attempt to close this gap by reviewing the foundations of DFT for open-shell 
systems. In our presentation, we will pay particular attention to the role of the spin- 
density in DFT and to possibilities for targeting different spin states within the exact 
theory. Even though we will not discuss the currently available approximations in detail, 
we believe that for the future development of better approximations, it is crucial to know 
which exact theory is to be approximated. This is also a prerequisite for deriving ex- 
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act conditions on the approximate functionals, for setting up model systems that can be 
treated exactly, and for obtaining benchmark results from accurate wave-function theory 
calculations. 



This work is organized as follows. First, we introduce spin in the context of nonrelativistic 
quantum chemistry in Section [2] This is followed by a discussion of Hohenberg-Kohn (HK) 
DFT, highlighting the role of the spin density and of spin states for open-shell systems 
in Section |3} Next, the treatment of spin in the Kohn-Sham (KS) framework of DFT 
is reviewed in Section |4} It turns out that different options exist for deriving KS-DFT 
for open-shell systems, which are discussed and compared in detail. For completeness, in 
Sec. [5] we discuss DFT within the relativistic framework where spin is no longer a good 
quantum number. Finally, some possible future directions for DFT applied to open-shell 
systems are outlined in Section |6| 



2 Spin in Nonrelativistic Quantum Chemistry 

2.1 Spin Structure of the One-Electron Wavefunction 

The nonrelativistic quantum-mechanical equation of motion for a single electron in an 
external electrostatic potential fext(^) is provided by the time-dependent Schrodinger 
equation (SE), which in Gaussian units reads, 

h^{r,t) = [f + q,v,,t{r)]^{r,t) = ih^^^{r,t), (1) 

with the kinetic energy operator T = = A, where p = — i/iV is the momen- 

tum operator and rrif. and qe are the mass and the charge of the electron, respectively. 
Stationary states can then be obtained from the time-independent Schrodinger equation. 
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(2) 



In a nonrelativistic framework, spin is introduced in an ad hoc fashion by employing a 

, i.e., 



two-component representation for the wavefunction 54[ 55 



(3) 



For a rigorous introduction of spin in quantum chemistry, it is necessary to start from 
relativistic quantum mechanics, where spin is naturally included in the Dirac equation. 
This will be discussed later on in Section [5l 

In the nonrelativistic two-component picture, spin-independent operators — such as the 
one-electron Hamiltonian h in Eqs. ([T]) and ([2]) — act on both of these components, i.e., 
they are proportional to the 2x2 unit matrix I2. The two-component structure of the 
wavefunction is only probed by operators expressed in terms of the Pauli matrices. 



I, 1^ 



(4) 



In particular, the operator corresponding to the electron spin is 
s = -a- = -[a^,ay,(x,) 



(5) 



The three components of this spin operator fulfill the same commutation relations as those 
of the angular momentum operator 1, i.e.. 



(6) 



which is the basis for considering spin as an intrinsic angular momentum vector s. For 
the squared magnitude of the electron spin, one obtains the diagonal operator 



"2 -2,-2,-2 ^ fe2 1 

S = S^ + Sy + S^ = -h h 



(7) 



and as for the angular momentum, this operator commutes with each component of the 
spin (i.e., [s^, Sa] = for a = x, y, z). 
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In nonrelativistic quantum mechanics, one postulates that the spin operator s is related 



to an intrinsic magnetic moment of the electron 54 , 56 , 57 . This spin magnetic moment 
is described by the operator 



h 



Us = = -hb cr, 



I I ^ 

where fin = — - — is the Bohr magneton, c is the speed of light in vacuum, and the factor 

znieC 

two is the electron (yf-factor when neglecting quantum electrodynamical effects. 

Instead of explicitly writing two-component wavefunctions and 2x2 matrix operators, in 
quantum chemistry it is common to employ a different notation, which will turn out to 
be particularly convenient for handling many-electron systems. Namely, as a shorthand 
notation, one introduces the orthonormal spin functions a{s) and which depend on 



a spin variable s 54,56,58 . This spin variable can only assume the values +| and — |, 
and the spin functions are defined such that 

a{+l) = 1 and = 0, (9) 

/3(+|) = and /3(-|) = l. (10) 

Then, the wavefunction of Eq. (jsj) can be expressed as 

ij{r,s)=MrMs)+Mr)m, (11) 

with the first component given by ip{r, +1) = ipa{f) and the second component given by 
ip{r, — ^) = ipp{r). It is important to realize that the spin functions a and (3 are merely 
a way of expressing two-component wavefunctions, in which the spin variable s has the 
role of labeling the different components. Spin-independent operators are then given by 
a one-component operator acting only on the parts of the wavefunction that depend on 
the spatial coordinate r, while the spin operators s^, Sy, s^, and act only on the parts 
depending on the spin variable s. 

The nonrelativistic one-electron Hamiltonian h is spin- independent and hence commutes 
with all spin operators, in particular [h, Sz] = [h, s^] = 0. Therefore, its eigenfunctions 



can be chosen as eigenfunctions of and of s^- The eigenfunctions of Sz are given by 
a{s) and (5{s) and the wavefunctions of an electron which are also eigenfunctions of Sz 
are thus of the form 

il){r , s) = il){r)a{s) or il){r,s) = il){r)(5{s). (12) 

For any spatial part of the wavefunction, these two different total wavefunctions are 
possible, i.e., each eigenvalue of the nonrelativistic one-electron Hamiltonian is two- fold 
degenerate. The first is identified with an a- or "spin-up" electron {sz = +h/2), whereas 
the second one corresponds to a /3- or "spin-down" electron {sz = —h/2). Of course, also 
linear combinations of these two degenerate eigenfunctions are solutions of the Schrodinger 
equation. Nevertheless, for a single electron any eigenfunction can be expressed as a 
product of a spatial part and a spin part. 



2.2 Spin Structure of the Many-Electron Wavefunction 

Within the two-component picture introduced in the previous section, the one-electron 
Hilbert space ^"H is spanned by all admissible one-electron wavefunctions (which have 
two components related to the spin of the electron). For an A^-electron system, the 
wavefunction is an element of the corresponding A^-electron Hilbert space ^"H, which is 
the tensor product space of the Hilbert spaces of each electron, i.e.. 



N 



n = ^n{i)®^ni2)®---^^n{N), (13) 

where the number given in parentheses designates the corresponding electron. As the 
one-particle wavefunctions each have two components, the A^-electron wavefunctions have 
2^-components, and operators have the dimension 2^ x 2^. For a more detailed discussion 



of the tensor structure of the many-electron wavefunction, see, e.g., Refs. 19-21 and 



chapter 8.4 in Ref. 57 



As an alternative to explicitly handling many-component wavefunctions, it is again con- 



venient to introduce spin coordinates Sj that can be used to distinguish the different 



components 56 , 59 , 60 . Then the A^-electron wavefunction depends on N spin coordi- 



nates in addition to the spatial coordinates, 

^ = ^(n, si, ra, S2, . . . , r-AT, sn) = ^(a^i, X2, . . . , xn), (14) 

where Xi = (rj, Sj) denotes the combination of spatial and spin coordinates. Each spin 
coordinate can assume the values —1/2 and +1/2. The possible combination of values 
for these spin variables each corresponds to one component of the many-electron wave- 
function. In total, 2^ different combinations of values are possible, that are thus used to 
label the 2^ components of the many-electron wavefunction. 

The nonrelativistic Hamiltonian describing N electrons in an external electrostatic po- 
tential fext(^) is (in Gaussian units) given by. 



N r ^0 N N n 

Te 



H = Y: -^A. + q^v^ir^ + E E T' ^'^^ 

1=1 -' 1=1 J=l+1 ■' 

where the Laplace operator Aj acts on the coordinate of the ith electron and rjj = |rj — rj| 
is the distance between electrons i and j. This Hamiltonian does not contain any spin- 
dependent terms and only acts on the spatial coordinates. 



The operator of the total spin S of a many-electron system is obtained 56 ,58 by summing 



the spins of the individual electrons S = X]iLi^('5i)- particular, one has for the z- 
component of the total spin operator, 

N 

Sz = ^Sz{si), (16) 

i=l 

and for the square of the total spin 

N N N N 

'^' = EE^(^')-^(^^) = i^^'+2E E ^(^o-^(^.), (17) 

i=l j=l i=l j=i+l 

where the first term emerges because any many-electron wavefunction is an eigenfunction 
of s^{si) with eigenvalue (3/4) (i.e., electrons are spin-1/2 particles). Note that S'^ 
couples different electrons, i.e., it is a two-electron operator [611. 



Both S"^ and 5*^ commute with the Hamiltonian of Eq. (15) and with each other, i.e., 

[H,S^] = [H,S,] = [S\S,]=0. (18) 

Therefore, the eigenfunctions \1/ of the Hamiltonian can always be chosen as eigenfunctions 
of Sz and with, 

5^^ = 5(5+1)^2^ (^g) 

Sz'^ = Msh^ with Ms = -S,...,+S. (20) 

In general, eigenfunctions of the Hamiltonian belonging to different eigenvalues of S"^ 
have different energies, while for each energy eigenvalue there are always 25+1 degen- 
erate eigenfunctions differing in M5. Since all three components of S commute with the 
Hamiltonian, any choice of the quantization axis is possible and will lead to identical 
results. 



Finally, the spin structure of the many-electron wavefunction 59 is also determined by 



the Pauli principle 62,63 . It requires that the wavefunction is antisymmetric (i.e., it 
has to change sign upon exchange of two electrons). This can be expressed with the 
permutation operator, 

"Pij^i- • • ) ^i; "Sj, . . . , Sj, . . . ) \E'(. . . , T^j, Sj, . . . , Sj, . . . ), (21) 

as Vij"^ = — \E'. Here, Vij exchanges both the spatial and the spin coordinates of elec- 
trons i and j. To express this requirement in a different form, one can introduce the 
antisymmetrizer A defined as 
m 



'^=47tE(-1)'^^' (22) 
p=i 

where the permutation operators Vp are ordered such that even numbers p are assigned 
to those that are generated by an even number of pair permutations, and odd numbers 
denote those generated by an odd number of pair permutations. Then, the requirement 
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that the wavefunction is antisymmetric with respect to any pair permutation is equivalent 
to requiring = -\/iVT\l/, i.e., ^ has to be an eigenfunction of A. 

Since this Hamiltonian does not contain terms that couple the spatial and the spin coor- 
dinates, one could naively expect that — as in the one-electron case — the many-electron 
wavefunction can always be written as a product of a part depending on the spatial coor- 
dinates and of a part depending on the spin coordinates. However, for systems with more 
than two electrons this is in general not the case. The antisymmetrizer A then contains 
a sum of permutation operators each acting on both spatial and spin coordinates. Hence, 
it couples spin and spatial coordinates so that its eigenfunctions cannot be expressed as 
a product of a spatial and a spin part. Thus, the structure of the many-electron wave- 
function with respect to the exchange of spatial coordinates is dependent on the spin 



structure 64 66 . This is the most important consequence of the presence of spin in a 



nonrelativistic theory. 

2.3 Spin Structure of the Electron Density and Spin Density 

The (total) electron density p{r) describes the probability density for finding any electron 
of a many-electron system at position r. It can be calculated from the wavefunction as 

p(r) = iV j \'^{r,si,X2,. . . ,xiy)\'^ dsidx2 - ■ -dxN, (23) 

i.e., by integrating the squared absolute value of the wavefunction over all but one spatial 
coordinate. By writing out the integration over the corresponding spin variable explicitly, 

p(r)=N j \'^{r,+l,r2,S2,. . . ,rN,SN)\^ d^r2ds2 - ■ -d^rNdsN 

+ N j |*(r, -I, ra, S2, . . . , r^, s^)P dVsdsa ■ ■ ■ dVjyds^ 

= p^{r)+p^{r), (24) 
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one notices that it is a sum of components Pa{f) and p/^ir) that can be interpreted as the 



probabihty densities of finding an a- or a /3-spin electron [56| 60 . Their integrals, 

N^ = jpair)d^r and Np = jpp{r)d\, (25) 
give the number of a- and /3-electrons. 

It is then natural to define the spin density, which gives the excess of a-electrons at a 
given point, as 

Q[r) = N j ^*(r, si, 032, • • • , xn) cx^(si) ^'(r, si, a^a, • • • , xn) dsida^a ■ ■ ■ da^TV 

= p^{r)-pp{r). (26) 

Here the Pauli matrix ^^(si) operates on the first spin coordinate only. Just as the 
electron density, which is probed by X-ray diffraction experiments, the spin density is 
an observable. The spin density Q{Ri) at the position of a nucleus / is accessible from 
electron paramagnetic resonance (EPR) experiments, where it determines the nuclear 



hyperfine coupling constants 56,67,68 . Similarly, the spin density at a nucleus can give 



rise to shifts in paramagnetic nuclear magnetic resonance (pNMR) l68 71 . Full spatially 



resolved spin densities can be determined in neutron scattering experiments 72 77 



The expectation value of a multiplicative one-electron operator such as V = Yl!i=i "^(^i) 
can be calculated directly from the electron density, 

ly) = {^\V\^) = j p{r)v{r)d^r, (27) 

i.e., the full wavefunction is not needed. Similarly, expectation values of spin-dependent 
operators expressed only in terms of cr^ can be obtained from the spin density. In partic- 
ular, the expectation value of is given by 

{S.) = ^ [ Qir)d\ (28) 



2 _ 

and for eigenfunctions of Sz, one has 

Ms = lJ Qir) dV = ^(iV. - N^). (29) 
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The contribution of the electron spin magnetic moments to the interaction of a molecule 
with an inhomogeneous external magnetic field Bg^t^r), the so-called spin Zeeman inter- 



action, is determined by the operator 56 



N N 



Hz = - ^ Be^tiri) ■ fisi-s-i) = fiB^ -Bext(n) " cr(si). (30) 



i=l i=l 



Hence, for an inhomogeneous magnetic field in z-direction, i.e., Be:^t{r) = (O, 0, ^^(r))"^, 
the expectation value of the spin Zeeman interaction can be evaluated directly from the 
spin density as 

{Hz)=fiB ! Q(r)5,(r)dV. (31) 



If one considers an eigenfunction of S"^ with eigenvalue S{S + 1)^^, this eigenvalue is 
(25* + l)-fold degenerate and one can construct a set of 25 + 1 eigenstates of with 
eigenvalues M^/i, where Ms = —S, . . . , +S. The total electron densities p^^^ and the spin 



densities Q^''^{r) of these 5*^ eigenstates are related to each other 60,78,79 : All 25' + 1 
states share the same total electron density, 

pA/s(^) =pMs=5(^) (32) 

and the spin densities are given by 

gA/s(^)= WgM.=5(^) (33) 



. 5 ^ 

where p^^^^{r) and Q^'^^^^{r) are the total electron density and the spin density of the 
state with highest M5, respectively. Hence, the spin densities have the same functional 
form and are connected by a simple scaling. It immediately follows that Q'^^^{r) = 
—Q^'^'^^{r) and that the spin density vanishes for states with Ms = 0. 



3 Spin in Hohenberg— Kohn DFT 



Traditionally, quantum chemistry sets out to calculate approximations to the many- 
electron wavefunction \l/ of a molecule in its ground state by minimizing the energy ex- 
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pectation value with respect to under the constraint that \1/ represents a normahzed 
and antisymmetric A^-electron wavefunction, i.e., 



Ea = mm{m^\H\m^) with (^^|^^) = 1 and A^^ = ^N\m^, (34) 
where the nonrelativistic Hamiltonian H within the Born-Oppenheimer approximation 



was given in Eq. (15). In molecular systems, the external potential fext(^) is given by 
the Coulomb potential of the nuclei, i.e., fext(^) = 'i^nuc(^') = —1e Yli ^//k ~ where 
the sum runs over all nuclei with charges Zi at positions -R/. Thus, the nonrelativistic 
molecular Hamiltonian assumes the form, 

AT -2 N N 2 ^ 

H = Y1 -^^i + E E — + E ^--(^') =T + Vee + Vnuc (35) 



According to this structure of the Hamiltonian, the energy expectation value in the above 
minimization is usually split up as 

(^1^1^) = (^|f 1^) + (^iKel^) + (^iKucI^). (36) 

However, the wavefunction itself is not directly needed for calculating these expectation 
values. The evaluation of the first term, corresponding to the kinetic energy, only requires 
the one-electron reduced density matrix (1-RDM), whereas the second term describing 
the electron-electron interaction can be calculated from the diagonal two-electron reduced 
density (2-RDM) matrix. Finally, the electron-nuclear attraction energy can be evaluated 
directly from the multiplicative operator of the electron-nuclei Coulomb interaction and 
from the electron density only as. 



(^iKucI^) = qej Pir)vnncir) d'r. (37) 
This can be exploited by performing the minimization with respect to the 2-RDM di- 



rectly 80 , but difficulties arise in enforcing that the 2-RDM corresponds to an actual 



antisymmetric A^-electron wavefunction. 



Density-functional theory (DFT) (50 53 provides the theoretical framework for calcu- 



lating the energy expectation value directly from the electron density p{r) only. The 
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foundations of this exact theory will be outlined in the following, focussing on its appli- 
cation to open-shell systems. 



3.1 Hohenberg— Kohn Theorems 



The first Hohenberg-Kohn theorem 81 states that for each electron density p(r) that can 
be obtained from a ground-state wavefunction (such densities are called f-representable 
densities), the external potential fext(^) that yields this electron density as the ground 



state when employed in the Hamiltonian of Eq. ( 15 ) is unique up to a constant. Therefore, 
this potential is a functional of the electron density and since it completely determines the 
Hamiltonian, also the ground-state wavefunction — which can in turn be determined by 
solving the corresponding Schrodinger equation — is a functional of the electron density. 
Furthermore, all observables of the system, in particular the total energy in a given nuclear 
potential fnuc(^)5 can be obtained from this wavefunction. This connection between the 
electron density and the total energy is illustrated in Fig. [T^. Therefore, there exists an 
energy functional E[p\ that relates the electron density to the total energy. 



According to the second HK theorem [81], the ground state energy Eq of a system of 
electrons in a given nuclear potential Vmxc{'r) can be determined by minimizing the total 
energy functional 

Eq = mm E[p], (38) 
p 

with 

E[p] = QeJ p{r) t;„uc(r) dV + Fhk[p], (39) 

under the constraint that p(r) integrates to electrons. The electron density for which 
this minimum is achieved is the ground-state electron density po(^)- In this equation, 
the total energy functional has been split into a system-specific part (the first term). 
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depending on the nuclear potential, and a system-independent part (the second term), 
which is called the universal Hohenberg-Kohn functional Fhk[p]- 



Following the Levy constrained-search formulation of DFT 82 , 83 , this universal HK 
functional is given by 



FhkM =min(^ T + Ke (40) 

where T and Vee are the operators of the kinetic energy and of the electron-electron repul- 
sion energy, respectively. The minimization runs over all wavefunctions that yield the 
target electron density p. From these wavefunctions, the one with the lowest expectation 
value of T + Vee is chosen. The minimization of E[p] will lead to the exact ground-state 
electron density po, and the exact ground-state wavefunction \&o is the one for which the 



minimum in Eq. (40) is obtained. This (nonrelativistic) ground-state wavefunction \I'o 
has to be an eigenfunction of S"^. Therefore, it is sufficient to restrict the constrained 
search to wavefunctions that are eigenf unctions of S^. 

Using the Levy constrained search for defining the HK functional also extends the domain 
in which the above functionals are defined from f-representable densities (i.e., densities 
that are obtained from a ground state wavefunction) to A^-representable densities (i.e., 
densities that are obtained from any wavefunction, not necessarily a ground state). The 
resulting generalization of the first HK theorem is illustrated in Fig. [TJo. More details 
on the foundations of HK-DFT and on the more general definition of the universal HK 
functional introduced by Lieb [84 1 can be found in dedicated reviews on these topics 
|5Tl[85l[86l. 



Initially, Hohenberg and Kohn explicitly excluded degeneracies in their derivation of the 



HK theorems 81 . However, a state with 5 > will be (25* + l)-fold degenerate, with 



the different degenerate wavefunctions \E'^^s corresponding to Ms = —S, . . . , S. The nec- 



essary generalization of the HK theorems is possible in a straightforward way 87 . The 



degenerate wavefunctions vE^^^s and all their linear combinations share the same electron 
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density. For a given p, the minimum in Eq. (40) is achieved for all degenerate wavefunc- 
tions spanned by the S^-eigenfunctions This is illustrated in Figure [l]^. Therefore, 

the minimization of the total energy functional E[p\ will still lead to a unique ground-state 
density po- If it is necessary to obtain a unique minimizing wavefunction, the constrained 
search can be restricted to wavefunctions corresponding to a specific value of Mg- In this 
case, the wavefunction is again uniquely determined by the electron density. 



3.2 Spin Density in Hohenberg— Kohn DFT 

According to the HK theorem only the total electron density is required for obtaining the 



exact ground-state energy and electron density 81,88 . Therefore, irrespective of the spin 
state, the spin density Q{r) or the individual a-electron and /3-electron densities Paij") 
and pp{r) are not required during the minimization of the total energy functional, and 
the ground-state spin density Qoir) is not directly available. However, if a suitable value 
of Ms is chosen, the wavefunction vlf^^s jg uniquely determined by the total density p{r), 
i.e., \Ef*^'S = \Ef^s|^pj_ From this wavefunction, the ground-state a-electron and /3-electron 
densities, 



[p\{r) = N J |^^^«[p](r,+i,r2,S2,...)|' d^dss ■ ■ ■ d Vds^ (41) 
P0'[p]{r) = N j |^*^«[p](r,-i,r2,S2,...)|' dV2ds2 ■ ■ ■ dV^vds^, (42) 

as well as the corresponding spin density 

Q'''[p\{r)=p'J^[p\{r)-pf^[p\{r) (43) 

can be calculated. In these expressions, the superscript Ms indicates that a specific 
value of Ms has to be selected if one is interested in individual a- and /3-densities or 
in the spin density. The chosen value of Ms fixes the number of a- and /3-electrons by 
Ms = 1{N''-N'^). 
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In order to generalize HK-DFT to use the individual a- and /3-densities instead of the 



total density only, the minimization of the total energy can be rewritten as 50[ 89 

mm -Iqe I p(?')fnuc(^) d^r + min 



Eo = min ii^[p] = min ■{ Qe pMvnucir) d r + min ( ^ 
p p \ J \ 











min ( ^! 




•>] 


) 









and by removing the outer minimization with respect to p, one obtains 



Eo = min <^ / [p'^{r) + p^(r)l v^ndr) dV + min 



(44) 



(45) 



This defines a universal HK functional in terms of the a-electron and /3-electron densities, 



i^HK[p",p1= min 



VP . 



(46) 



In contrast to the HK functional of Eq. (40), the constrained search in this spin- resolved 



HK functional now runs over all wavefunctions corresponding to a given pair of p" and 
. By minimizing the generalized total energy functional [cf. Eq. (45)], 



E[p",p^] = qej [p"{r)+p^{r)]vnUr) dV + Fhk[p", P^], 



(47) 



under the constraint that p° and p^ integrate to A^" and A^'^ electrons, respectively, it is 
then possible to obtain the ground-state a-electron and /3-electron densities pg and pg. 
Again, in this minimization a specific value of Ms = ^(A^" — A^^) has to be selected. 

Instead of using p" and p^, it is also possible to employ the total electron density p = 
p" + p'^ and the spin-density Q = p" — p^ as variables of the spin-resolved HK functional. 
This gives 



^hk[p,Q] = min 



(4J 



While usually — in particular in practical applications of spin-DFT — it is more common 
to employ the a- and /3-densities as basic variables, in the following such a formulation in 
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terms of p and Q will often turn out to be useful, because it allows for an easier comparison 
to the spin-independent functionals defined only in terms of the density p. We will switch 
between these two representations whenever suitable. 

With the spin- resolved HK functional Fhk[p, Q], it is now also possible to give a simpler 
prescription for obtaining the spin-density corresponding to a given total density: Q*^s[p] 
is the spin density for which Fhk[p, Q] is minimized, under the constraint that Q^^ 
integrates to twice the chosen value of M5, i.e., 

Q"'^ [p] = arg min Fhk [p, Q'''] with I [ Q''' (r) dV = M5. (49) 
qMs 2 J 

Therefore, the spin- independent HK functional F[p] can be obtained from the spin- 
dependent HK functional Fhk [p, Q] as 

Fhk[p] = Fhk[p, Q[p]] = min Fhk[p, Q], (50) 

Q 

where Q[p] is any of the spin densities that minimize -Fhk[P) Q] for the given total density 
p. Of course, as long as no specific M5 is chosen, Q[p] is not unique, but any admissible 
choice must lead to the same energy. 

Finally, we have to consider whether there exist extensions of the HK theorems that justify 
the use of p and Q (or of p" and p'^) as basic variables. For a generalization of this kind, 
it is not sufficient to consider wavefunctions generated by an external potential Vext{r), 
but also wavefunctions obtained in the presence of an additional external magnetic field 
Byir) have to be taken into account. Such an extension of the HK theorems was first 



given by von Barth and Hedin 90 , and was only recently put on more firm ground by 



extensions of Lieb's formulation of DFT 91 92 . Similar to the external potential, which 



is only known up to a constant, it is also possible to add a constant shift to the external 



magnetic field Bz{r) without changing the wavefunction or the (spin-) density [93 94 



This leads to a number of peculiarities related to the differentiability of the spin-dependent 



energy functional 95 ,96]. However, most of these issues do not appear if the treatment 
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is restricted to eigenfunctions of (the case of interest here) [97] or can be addressed by- 



constraining Ms to a fixed value (as we are always requiring here) [98, 99 



3.3 Fractional Spin Conditions on the HK Functional 

For states that are not a singlet (i.e., for S > 0), there are different degenerate wavefunc- 
tions \E'*^'5 [p] . These wavefunctions all share the same electron density, but correspond to 



different spin densities Q'^^^[p]- These different spin-densities are related by [cf. Eq. (33) 



Q^^[p](-)=('^)Q^^-^[p](r), (51) 



S 

and it follows that Q^^^{r) = -Q~^^{r) and that one obtains Q^^s=0{r) = for even 
values of S. For all these degenerate spin-densities, the exact total energy functional 
E[p,Q] must yield the same value. 

This statement can be generalized to linear combinations of these spin-densities. Within 



an ensemble formulation of spin-DFT, it can be shown 100 that for any properly normal- 
ized linear combination of the spin-densities Q^^^ [p] , the same energy should be obtained. 
Therefore, one finds that 



E[p,aQ^'^^=^[p]] = const, for - 1 < a < +1, (52) 

which implies for the spin-dependent HK functional 

Fhk [p, a Q*^^=^[p]] = const, for - 1 < a < +1. (53) 

Such ensemble spin densities correspond to a situation with a non-integer number of a- 
and /3-electrons {^^ fractional spins"). As has been pointed out by Yang and co-workers, the 
above constancy conditions for fractional spins, which is a property of the exact energy 
functional, is violated by all contemporary approximations. Therefore, it was suggested 
that many problems appearing in practical DFT calculations with such approximations 



might be connected to this violation 49, 101 102 
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The use of fractional spins makes it possible to further simplify the relation between 



the spin-dependent and spin- independent HK functional given in Eq. (50), because one 
realizes that 

i^HK[p] = i^HK[p,Q = 0]. (54) 

This equation holds both for systems with an even number of electrons (where one can 
always choose Ms = 0, corresponding to Q{r) = 0) and for systems with an odd number 
of electrons (where Q{r) = is only possible if one allows for fractional spins). 



3.4 Spin States in Hohenberg-Kohn DFT 



In their initial formulation 81 , the HK theorems were applicable only to the ground state. 
In particular, minimizing the total energy functional E[p] only yields the ground-state 
energy and electron density (and if the spin-dependent energy functional is employed, also 
the spin density). However, as was first shown by Gunnarson and Lundquist, HK-DFT 



can be generalized to the lowest-energy states of a given symmetry 103 . Of particular 
importance is the calculation of the lowest state of a particular spin symmetry, i.e., of 
the lowest state with a particular eigenvalue of S'^. Such a generalization is most easily 
presented within the constrained-search formulation of DFT. 

In order to obtain the energy and electron density of the lowest state corresponding to a 
given value of S, one has to define the spin-state specific energy functional, 

E'[p] = qej p{r) i;nuc(r') dV + F|k[p], (55) 

with the spin-state specific HK functional, 

F^kH = min (^^ f + Ke with S^^^ = S{S + l)ff^^, (56) 

where the constrained search now only includes wavefunctions \I''^ which are eigenfunctions 
of 52 with the proper eigenvalue. Therefore, one obtains a different HK functional and 
thus a different total energy functional for each value of S. 

21 



Within spin-DFT, the simplest way of obtaining a functional that at least partly allows 
one to select certain spin states is by choosing an appropriate value for Ms- This way, 
only states with S > Ms are accessible because for spin states corresponding to a smaller 
values of S, the chosen value of Ms is not admissible. By minimizing the energy functional 
Ev[p,Q] under the constraint that Q{r) integrates to 2Ms, the lowest-energy state with 
S > Ms is obtained, i.e., 

E^^*-^-[p]=E[p,Q^-^-[p]]. (57) 

Equivalently, one can, of course, also minimize Ey[p°',p^] under appropriate constraints 
for A^" and A^'^. However, fixing Ms to target a specific spin-state is not completely 
general since the minimization is only restricted to states with S > Ms, not to states 
with a specific S. While it is, for instance, possible to calculate the lowest triplet {S = 1) 
state if the ground-state is a singlet {S = 0), it is not possible to target the lowest singlet 
state if the ground-state is a triplet. Therefore, to be able to calculate the lowest state of 
a given S, it would in general be necessary to employ the true spin-state specific energy 



functional of Eq. (55). 



4 Spin in Kohn-Sham DFT 

While the Hohenberg-Kohn formulation of DFT is exact, it is very difficult to set up 
computationally feasible, but nevertheless accurate approximate realizations of it. This 
is mainly rooted in the difficulty of approximating the kinetic-energy contribution to the 



HK functional as a functional of the electron density only 104 , 105 



A possible way out of this dilemma, that forms the basis of almost every present-day 



application of (approximate) DFT calculations, was suggested by Kohn and Sham 106 
Instead of considering the kinetic energy of the true system of interacting electrons, they 
proposed to calculate the kinetic energy of a reference system of noninteracting electrons 
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with the same electron density instead. This then aheady accounts for the largest part of 
the kinetic energy, and only a small remainder has to be approximated. The KS approach 
still allows for the formulation of an exact theory, which will be outlined in this section. 

In KS-DFT one considers two different quantum-mechanical systems at the same time: 
The true molecular system of interacting electrons and a reference system of noninter- 
acting electrons. The link between these two systems is established by requiring that 
their electron densities p{r) and Ps{i") are equal (see Fig. |2|. Their wavefunctions, how- 
ever, will in general be different. For open-shell molecules, different options exist for 
introducing such a reference system: The first option is to require only that the electron 
densities of the interacting and the noninteracting systems agree. This leads to a spin- 
restricted KS-DFT formulation. The second option is to require that in addition to the 
total electron densities, also the spin densities of the two systems agree. This results in a 
spin-unrestricted formulation of KS-DFT. Of course, these two options are equivalent for 
closed-shell systems (i.e., for singlet states with 3 = 0). 

Note that any version of KS-DFT relies on the assumption that such a noninteracting 
reference system with the same electron density (and possibly also the same spin density) 
as the interacting system exists. In practice, it is always assumed that this so-called 
fs-representability condition is fulfilled, even though this is not guaranteed and several 



counter-examples are known 84 ,107 111 . For a detailed discussion of these subtle issues. 



see, e.g., Refs. 85 . 
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4.1 Spin-Restricted Kohn-Sham DFT 



4.1.1 Noninteracting Reference System 

First, we consider a system of N noninteracting electrons in an external potential Vs{r). 
Such a system is described by the Hamiltonian 

^ fc2 ^ 

1=1 1=1 

where the subscript s (for "single-particle") is introduced to indicate that the quanti- 
ties refer to a system of noninteracting electrons. Because this Hamiltonian does not 
contain terms that couple different electrons, an exact wavefunction is given by an anti- 
symmetrized product of one-electron functions (i.e., by a single Slater determinant), for 
which the short-hand notation 

^^(cci, . . . ,ccAr) = (pia,^pil3,(p2a,^p2P,. . . (59) 

can be used. The spatial parts of the one-electron functions (orbitals) can be obtained as 
the solutions of the one-electron equation 



2m, 



(Pi{r) = ei(pi{r), (60) 



which trivially emerge from the energy eigenvalue equation Hg^s = Eg^s with Eg = e^. 
Since this one-electron Hamiltonian does not depend on the spin of the electron, each 
energy eigenvalue is two- fold degenerate (i.e., each spatial orbital can be combined with 
an a- or with a /3-spin function). In particular, the spatial orbitals are identical for a- 
and for /3-electrons, and the resulting Slater determinant is therefore spin-restricted. 

In such a spin-restricted Slater determinant, each spatial orbital (pi can either be doubly 
occupied (i.e., it appears both in combination with an a- and with a ^-spin function) or it 
can be singly occupied with either an a- or a /3-electron. For the ground state, such singly 
occupied orbitals can only occur for the highest occupied molecular orbital (HOMO), and 
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more than one singly occupied orbital can only be present if the HOMO is degenerate (for 



a detailed discussion, see, e.g., Ref. 109 and chapters 3.3 and 3.4 in Ref. |53|). 
For a single, spin-restricted Slater determinant, the electron density is given by 

occ. 

p{r) = J2f.\v.{r)\^ (61) 

i 

where the occupation numbers fi are either 1 or 2. The spin density is determined only 
by the singly occupied orbitals, and can be calculated as 

singly occ. 

Qsir)= ^^Iv'^Wr, (62) 

i 

where the sum only runs over the singly occupied orbitals and where Sj = +1 for a-spin 
and Si = — 1 for /3-spin orbitals, as all doubly-occupied orbital contributions drop out 
because of the identical spatial distribution of these a, /3-pairings. Finally, the kinetic 
energy of a spin-restricted Slater determinant can be calculated as 

occ. occ. « 

Ts = Y,f^{v^\f\v^) = -^Y.f^^ J V>^ir)A^^ir) d'v. (63) 

Thus, it turns out that spin-restricted Slater determinants in which the occupation num- 
bers fi are identical share the same electron density and have the same kinetic energy. 
Such determinants — which can only differ in the spin of the singly occupied orbitals — 



are degenerate with respect to the noninteracting Hamiltonian of Eq. (58) 



The noninteracting Hamiltonian Hg commutes with both Sz and with S"^. Therefore, it 
is always possible to combine degenerate eigenf unctions such that they are also eigen- 
functions of Sz and S"^. A restricted Slater determinant is always an eigenfunction of Sz 



with eigenvalue Mgh = {N^ — Np)h/2 [58,60 . If all singly occupied orbitals are either 
a- or /3-spin orbitals, then it is also an eigenfunction of S"^ with S = \Ms\. In all other 
cases, an eigenfunction of S"^ can be constructed as a linear combination of (degenerate) 
determinants in which the same orbitals are singly occupied and which correspond to the 
same value of Ms- Such linear combinations are known as configuration state functions 
(CSF) [To]. It is important to understand that for any eigenfunction of the noninteracting 
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Hamiltonian that is an eigenfunction of S"^ with eigenvalue 5'max(5'max + 1)^^, degener- 
ate eigenfunctions with S = 0, 1, ... , 5'max (for an even number of electrons) or with 
S = ^, |, . . . , Sraax {^OT an odd number of electrons) can also be constructed. Thus, for 
each energy eigenvalue, there is one CSF corresponding to S = (i.e., a singlet state) or 
S = ^ (i.e., a doublet state) for an even or odd number of electrons, respectively. 

The Hohenberg-Kohn theorems still hold for a system of noninteracting electrons. Thus, 
the ground-state density of Hg can be determined by minimizing the noninteracting energy 
functional, 



E,[p] = T,[p] + qj pir)vsir) dV, 



(64) 



where the noninteracting kinetic-energy functional Ts[p] can be defined in the Levy con- 
strained-search formalism as 

Ts[p] = mm {^s\f\^s). (65) 

In this definition, the constrained search includes all wavefunctions \E's that correspond to 
a system of noninteracting electrons with density p. As discussed above, this could be fur- 
ther restricted to singlet or doublet wavefunctions. In this definition, [p] is independent 
of the spin density. 

The ground-state electron density po is obtained from minimizing Es[p] under the con- 
straint that the number of electrons is preserved, and the corresponding ground-state 



wavefunction \l/s,o is the one for which the minimum in Eq. (65) is achieved. Again, this 
wavefunction could always be chosen as a singlet or a doublet for an even or odd number 
of electrons, respectively. 
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4.1.2 Interacting Energy Functional and Exchange— Correlation Energy 



The (spin-resolved) HK functional of the true system of interacting electrons [Eq. (48) 
can now be decomposed as 

FnK [p, Q] = Ts [p] + J[p]+ E,, [p, Q] , (66) 

where Ts[p\ is the noninteracting kinetic energy introduced in the previous section, J[p] 
is the classical Coulomb interaction of the electron density with itself, 

^W = f/4^dVdV. (67) 
2 J \r — r'\ 

and the exchange-correlation energy E^dp^Q] is defined to account for the remaining 
energy contributions 

Exc [p, Q] = FuK [p, Q] - Ts [p] - J[p] . (68) 

This exchange-correlation functional could also be expressed in terms of the a- and /5- 
electron densities. The noninteracting kinetic energy Ts[p] is different from the true kinetic 
energy of the fully interacting system T[p]. Therefore, the exchange-correlation energy 
also contains the difference Tc[p] = T[p] — Ts[p] between the kinetic energy of the true 
interacting system and the kinetic energy of the noninteracting reference system. 

With these definitions, the total energy functional of the interacting system can be ex- 
pressed as 



E[p, Q] = T,[p] + J[p] + E,,[p, Q] + p(r)i;nuc(r')dV. 



(69) 



The spin-independent analogues of the exchange-correlation and the total energy func- 
tional, E^c[p] and E[p], are recovered from these definitions when setting Q{r) = [cf. 



Eq. (54);. 



The ground-state density of the true interacting system can be determined by minimizing 
the total energy functional E[p] with respect to p, under the constraint that it integrates 
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to the correct number of electrons. With the exact functionals, this will lead to the exact 
ground-state electron density Po{r). The corresponding exact wavefunction \E'o is the one 



for which the minimum in the constraint search in Eq. (48) is obtained. This ground-state 
wavefunction has to be an eigenfunction of and Sz, and the corresponding value of 5* 
determines the spin multiplicity of the ground-state, whereas the different Mg-states are 
degenerate. In contrast to the noninteracting case, the ground-state wavefunction is not 
necessarily a singlet or a doublet wavefunction, but could have a higher spin multiplicity. 

4.1.3 Kohn— Sham Potential 



Minimization of the noninteracting energy functional [Eq. (64)] with respect to the total 
density p, under the constraint that p integrates to electrons, yields the following 
Euler-Lagrange equation ISO], 



On the other hand, the total energy functional of the interacting system is given by 



Eq. (69), which upon minimization leads to the condition 



° " ^ ^ ^^^"""^^""^ + ^Coui[p]{r) + v,,[p]{r)) - P, (71) 

where fcoui[p](^) = Qe J p{f')/\f — f'\ d'^r' is the classical Coulomb potential of the elec- 
trons and v^c[p]{f) = i'^/le) SExc[p]/Sp{r) is the exchange-correlation potential. 

Since the definition of the noninteracting kinetic energy Ts[p] is the same in both min- 
imizations, and because we require that the total densities obtained for the interacting 
and the noninteracting system agree, we obtain for the KS potential 

^s[p]ir) = Vc^tir) + vconi[p]ir) + v^c[p]{r). (72) 

The ground-state electron density of the fully interacting system can thus be determined 
by solving the Schrodinger equation of a noninteracting system [i.e. with the Hamiltonian 
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of Eq. (58)) including the external potential Vs{r) given by Eq. (72)], and the orbitals of 
this noninteracting system can be obtained from the KS equations 



-A + qe{vnuc{r) + vcoui[p]{r) + v^c[p]{r)) 



ifi{r) = ei(pi{r). (73) 



2?7ie 

Thus, the ground-state electron density of the true interacting system is obtained from 
the wavefunction q of the noninteracting reference system. However, this ground-state 
wavefunction o of the noninteracting reference system does not agree with the ground- 
state wavefunction \l/o of the interacting system. Moreover, q can always be chosen as 
a singlet or doublet state (i.e., S = or S = |), whereas \l/o can correspond to any value 
of S. Therefore, the spin multiplicities of the ground-states of the noninteracting and of 
the interacting system can be different. 

4.1.4 Spin Density 

For S > 0, the ground state of the true interacting system is degenerate and there is 
a set of eigenfunctions of with different eigenvalues Mg. These wavefunctions ^I/^"^ 
have different spin-densities Q^'^lr) that are related by Eq. (33). By construction, the 



ground-state wavefunction q of the noninteracting reference system and the ground- 
state wavefunction \l/o of the fully interacting system only share the same electron density. 
However, the corresponding spin densities Q{r) and Qsi^) are in general not equal and 
the true ground-state spin density cannot be calculated from o- 



From Eq. (62) it is obvious that Q{r) and Qs(^) have to be different: For the spin- 
restricted noninteracting reference system, the spin density is determined only by the 
singly occupied orbitals and will thus have the same sign at every point in space (i.e., 
Qsi^) > for Ms > 0). However, it is known both from accurate calculations and from 
experiment, that for the interacting system the spin density has different signs in different 



regions in space 112 114 . 



To obtain the spin density in a restricted KS-DFT formulation, one has to minimize the 
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spin-resolved HK functional -Fhk[po) Q] defined in Eq. (48) with respect to Q(r), under the 
constraint that the spin density integrates to 2Ms. Since the exchange-correlation energy 
is the only part of this functional that depends on the spin density, the minimization with 
respect to Q{r) leads to the condition 

6Q{r) 6Q{r) " °' ^^^^ 



where the Lagrange multiplier A is zero because of Eq. (33). This suggests a two-step pro- 
cedure for determining the spin density in restricted KS-DFT. First, the total ground-state 
density po{f) is determined by solving the KS equations. Subsequently, the corresponding 
ground-state spin density Qo{r) can be calculated from the above minimization condition 
for a chosen value of Ms- 



4.2 Spin-Unrestricted Kohn-Sham DFT 



4.2.1 Noninteracting Reference System 



The choice of a reference system of noninteracting electrons with the same total electron 
density as the interacting system is not the only option. Alternatively, it is also possible 
to envisage a reference system of noninteracting electrons that has the same a-electron 



and /3-electron densities as the interacting system 90 ,103 ,112 . In this case, a reference 
system with the Hamiltonian 

i=l 1=1 



TV 



n 



N 



2m ~' ' -=A^L^"^' ""^^" ^.v ^/-v-»/j (75) 

i=l i=l 

is used. To distinguish them from those introduced earlier for a spin-restricted refer- 
ence system, the superscript "(u)" will be used for quantities referring to this spin- 
unrestricted reference system. Different potentials vf{r) = f*°*(r) + vf^'^{r) and ff(r') = 
vT^if) —vf^^^{r) for the a- and /3-electrons, respectively, are now needed in order to allow 
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the reference system of noninteracting electrons to have the same spin density as the in- 



teracting one. This corresponds to introducing an inhomogeneous external magnetic field 
in z-direction Bz{r) = —vf^"^{r) that only interacts with the electronic spins (i.e., the 



interaction due to orbital angular momentum is ignored) [cf. Eq. (30)]. 



An exact solution to the corresponding Schrodinger equation has the form of a single 
Slater determinant, but in contrast to the spin- restricted case the spatial orbitals now 
differ for a- and /3-electrons, i.e.. 



<l>(-")(a^i,...,x^) 



/3, v?2 



/3, 



(76) 



The spatial parts of the orbitals can be obtained from two separate sets of one-electron 
equations 



2m, 



2m, 



■A + ge<(r) 



and 



(77) 



Both the resulting a- and /3-orbitals form an orthonormal set, (v^^lv^") = Sij and {ip^\ipj) = 
6ij, but a- and /3-orbitals are in general not orthogonal to each other, i.e., (ipflipj) 7^ 0. 

The noninteracting Hamiltonian i^i"'* still commutes with S^, and any spin- unrestricted 
Slater determinant is an eigenfunction of Sz, with the eigenvalue Ms being determined 
by the number of a- and /3-electrons. However, in contrast to the spin-restricted case, the 
different eigenstates are not degenerate anymore. For constructing the ground-state 
wavefunction, the orbitals with the lowest orbital energies have to be occupied, which 
automatically fixes M5. By occupying other {a- or /3-electron) orbitals, excited state 
wavefunctions for the noninteracting reference system corresponding to different values of 
Ms can be obtained. 

However, ifi"^ does in general not commute with S'^, and the ground-state wavefunction 
is thus not an eigenfunction of anymore. Instead, the expectation value of can be 
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calculated as (assuming Na > [61,115 



JVc Np 



{S') = Ms{Ms + i)h' + h'Np-h'Y^Y. I 

i=l i=l J 



<(r)^?(r)dV 



(7J 



In the spin-restricted case, where the a- and /3-orbitals are equal and therefore mutually 
orthogonal, the last term equals the number of doubly occupied orbitals, and one obtains 
(S^) = Ms{Ms + l)h^. In the unrestricted case, this cancellation is only partial and a 
larger expectation value is obtained. This is often referred to as spin contamination. 



For an unrestricted Slater determinant the total electron density is given by. 



N„ 



i=l i=l 

and the spin density can be calculated as 



(79) 



(80) 



1=1 



1=1 



In contrast to the spin-restricted case, the spin density can now have different signs at 
different points in space. The kinetic energy of the unrestricted Slater determinant 
is 

.2 



Ti") = {^,\f\^,) =-^Y^I ^nr')A^r(r)dV 

^ i=l 



(81) 



and a noninteracting kinetic-energy functional can now be defined as 



min {^>^^^\f\¥^^) 



^2) 



or as 



Ti^\p,Q]= min (vl>W |f |vl>H). 



^3) 



In contrast to the spin-restricted case, this functional depends not only on the total 
electron density p(r) = Pai^) + Pis^r), but also on the spin density Q{r) = Pa{f) ^ Pp{'^)- 
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Usually, Ts^^ [p, Q] yields different kinetic energies for systems that share the same total 
electron density, but have different spin densities. 

In terms of a- and /3-electron densities, the total energy functional of the noninteracting 
system is given by 

Ei^\p^,p^] = T^^\pa,P,s] +qe[ Pair)vfir) d'r + qJ p^ir)v^ir) dV (84) 



The ground-state a- and /3-densities Po(^) Po(^) noninteracting system can 

then be determined by minimizing this energy functional with respect to pa and pp, under 
the constraint that these integrate to the correct number of a- and /3-electrons. 

4.2.2 Exchange— Correlation Energy Functional 



The spin-resolved HK functional (cf. Eq. (46)) of the true system of interacting electrons 
can now be decomposed as 

FnKW,P^] = T^'''\pa,pp] + J[p] + E^^[p^,p^], (85) 

where the spin-resolved exchange-correlation energy is defined as 

E^\Pa,Pp] = FnK[Pa,Pp] - Ti")[p„,p^] - J[p]. (86) 

Since this functional E^\pa, pp] has been defined via the spin-unrestricted reference sys- 
tem, it is in general different from the functional E^dPa, Pis] defined in the previous section. 
This difference arises because different definitions of the noninteracting kinetic energy are 
used in the two cases. 

With this definition of the exchange-correlation energy, the total energy functional of the 
true system of interacting electrons is given by 

E[pa, Pp] = Ti") pf,] + J[p] + Ei^^ [p^, pp] + q,j p{r)v^Ur) dV. (87) 
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This total energy functional is identical to the spin-resolved version of the total energy 



functional derived in the spin-restricted case [cf. Eq. (69)], even though it is decomposed 
in a different fashion. 

The ground-state a- and /3-electron densities p^iv^ and Po(^) then be determined by 
minimizing this total energy functional with respect to pa and p^, under the constraint 
that these integrate to A^^^ and electrons, respectively. Note that by choosing Na 
and A^^, a specific value of Ms is selected. As long as this Ms can be realized for the 
exact ground state, the resulting total ground-state density po(^) = Po (^) + Po(^) "^ill be 
independent of the choice of Ms- In addition, the minimization will then yield the exact 
ground-state spin density Q^^ir) = Po(r) — PqIv). 

4.2.3 Kohn-Sham Potential 

The minimization of the total energy functional Ei^\pa, p/s] of the spin-unrestricted non- 
interacting reference system leads to these Euler-Lagrange equations 

^^^^-^'"-^l^^ + ^^'CM-A^^O (88) 

dpa{r) dpa{r) 

For the interacting system, the minimization of E[pa,pp\ with respect to the a- and 
/3-electron densities yields 

^^}p2rf - = ^^^2l{r)''^ ^ ^^^''""^^''^ ^ ^Coui[p](r) + <K,P,](r)) - p„ = 0, 

(90) 

-^^, = ^4^^ + ^e(^nuc(r) + .coui[p](r) + Jp., P,] (r)) - p, = 0. 

(91) 

with the spin components of the exchange-correlation potential 

ar U ^ 1 5^ic^[Pa,P/3] , Rr . I 6E^^J[p^, pp] 

^xc Pa,P/3 r = and < Pa,P,3 r = ^ . 92 
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If we require that pa{f) and Pi3{r) — and thus both the total and the spin density of 
the ground state — are the same for the noninteracting reference system and the true 
interacting system, we find that the spin components of the KS potential are given by 

<[p^/](r) = vUr) + vcoM{r) + v'^^^ip^ , pP]{r) (93) 
vl[p'',P^]{r) = v,^,{r)+vconAp\{r) + v^^J,p'^,p^]{r). (94) 



Therefore, the exact ground-state a- and /3-electron densities of the true interacting sys- 
tem can be calculating by solving the Schrodinger equation of an auxilliary system of 



noninteracting electrons with the Hamiltonian of Eq. (75). The ground-state wavefunc- 
tion \E'^"J of this KS reference system is given by an unrestricted Slater determinant, 
constructed from the orbitals obtained from the KS equations. 



2m, 



-A + qe{vr,uc{r) + vc^^i[p\{r) + v'^^[p^, pp]{r)) 
■A + ge(^^nuc(r') + vcon\[p]{,r) + v^^[p^,pp\{r)) 



Vfir) = ef^f (r). (95) 



2me 

Here, the equations for the a- and /3-orbitals are coupled through the Coulomb and 
exchange-correlation potentials. 

Equivalently, the KS potential can be expressed as a component that acts on the total 
electron density, 

vT'lp, Q]{r) = ^ (<(r) + v^{r)) = v^r) + vcouM{r) + v'^^lp, Q]{r), (96) 
and one that acts on the spin density, 

vriP, Q]ir) = I (<(r) - v^ir)) = v^Ip, Q]ir), (97) 
where the total and spin exchange-correlation potential are given by 

These expressions will be used in the following section for comparing with the spin- 
restricted theory. 
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Even though the electron density and the spin density calculated from this J ^^^^ equal 
to those of the fully interacting system, it is important to realize that \l/^"o does not agree 
with the ground-state wavefunction \l/o of the true interacting system. In particular, \l/o 
can always be chosen as an eigenfunction of S^, whereas by construction, \1/^"q is not an 
eigenfunction of S"^ for S > 0. Thus, within an exact formulation of unrestricted KS-DFT, 
the wavefunction of the KS reference system is always spin contaminated for 5 > 0. 



4.3 Comparison of Restricted and Unrestricted Formulation 

The restricted and the unrestricted formulation of KS-DFT are based on different defini- 
tions of the noninteracting reference system for open-shell systems. In the spin-restricted 
case, the reference system is chosen such that its total electron density Ps{r) agrees with 
the one of the fully interacting system, while its spin density Qs{r) usually differs from 
the one of the interacting system. On the other hand, in the spin-unrestricted case the 
reference system is defined such that both its total electron density and its spin density 
agree with those of the fully interacting system. 

These different definitions of the noninteracting reference system have implications for 
the treatment of spin in KS-DFT. In the spin-restricted case, the wavefunction of 
the noninteracting reference system can always be chosen as an eigenfunction of S^. 
Nevertheless, the corresponding eigenvalue {S"^) = S{S + l)h'^ does not necessarily agree 
with the one obtained for the true interacting system. However, it is possible to require this 
equality with an additional constraint on the noninteracting reference system. In the spin- 
unrestricted case, the wavefunction of the reference system is not an eigenfunction of 
S"^, i.e., it is spin-contaminated. This is a direct consequence of the requirement that the 
correct spin density is obtained. Thus, the expectation value of S"^ becomes a complicated 



functional of the electron density 115,116 . Of course, the exact ground state density will 
still correspond to an interacting wavefunction that is an eigenfunction of S"^. 
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In both the restricted and in the unrestricted case, the wavefunction of the noninteract- 
ing reference system is an eigenf unction of Sz- Only in the spin-unrestricted case it is 
guaranteed that the corresponding eigenvalue Ms is the same as for the fully interacting 
system, but also in the spin-restricted case it can be chosen accordingly. These differ- 
ences between the restricted and the unrestricted formulation are summarized in Table HI 
One important observation is that it is impossible to set up a KS-DFT formalism such 
that for the noninteracting reference system one obtains both the correct spin density 
and a wavefunction that is an eigenfunction of (see also the discussion of this issue in 
Refs. (88|[TT2l). 



The different definitions of the noninteracting reference system in the restricted and the 
unrestricted formulations of KS-DFT also imply different definitions of the noninteracting 
kinetic energy, the exchange-correlation energy, and the exchange-correlation potential. 
These definitions are collected in Table |Tlj First of all, the use of different reference systems 
leads to different definitions of the noninteracting kinetic energy. In the spin-restricted 
case, Ts[p] is defined as the kinetic energy of a system of noninteracting electrons with 
the total electron density p{r) and is independent of the spin density Q{r). In contrast, 
in the spin-unrestricted case Ts^\p,Q] is defined as the kinetic energy of a system of 
noninteracting electrons with the total electron density p{r) and the spin density Q{r). 
These differ by the "unrestricted" contribution to the noninteracting kinetic energy, 

T^[p,Q] = T,[p]-T^-\p,Q]. (99) 

Only if the spin density vanishes, the restricted and the unrestricted definitions of the 
noninteracting kinetic energy are identical, i.e., Tu[p,Q = 0] = 0. 

Because of these different definitions of the noninteracting kinetic energy, a different 
decomposition of the HK functional is introduced in the restricted and unrestricted 
formalisms, respectively, which in turn leads to different definitions of the exchange- 
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correlation energy. These are related by 



(100) 



One important difference between the two formalisms is that in spin-restricted KS-DFT, 
the exchange-correlation energy E^dPiQ] is the only contribution to the HK functional 
that depends on the spin density, whereas in the spin-unrestricted theory both the exchange- 
correlation energy [p, Q] and the noninteracting kinetic energy Ti"^ [p, Q] depend on 



the spin density. Therefore, the fractional spin condition of Eq. (53), formulated for the 



HK functional in Sec. |3.3[ leads to different exact conditions for the exchange-correlation 
functional. In the spin-restricted case, the fractional spin condition applies directly to the 
exchange-correlation energy, 

^xc [p, a Q*^*'=^[p]] = const, for - 1 < a < +1, (101) 

whereas in the spin-unrestricted case it applies to the sum of the exchange-correlation 
energy and the noninteracting kinetic energy, 

Ti") [p, a Q'^'^=^[p]] + ^ie^ [p, a Q'"''^=^[p]] = const, for - 1 < a < +1. (102) 

Finally, the exchange-correlation potential (and thus also the resulting KS potential) dif- 
fers in the two formalisms. In spin-restricted KS-DFT, the exchange-correlation potential 
v^cIp] depends only on the total electron density and acts on electrons of both spin. On 
the other hand, in spin-unrestricted KS-DFT the exchange-correlation potential is differ- 
ent for a- and /3-electrons, i.e., it has two distinct components. The component of the 
exchange-correlation potential acting on the total electron density is given by, 
^tot^p g](^) = 1 SEt\p,Q] ^ 1 f 6E^,[p,Q] ^ ST^[p,Q] 



Qe 5p{r) Qe \ 5 p{r) 5p{r) 
= v^,[p,Q]{r) + Vu[p,Q]{r). (103) 

Here, the first term is the exchange-correlation potential in the spin-restricted formalism, 
while the second term Vu[p,Q]{r) = (l/qe) 6Tu[p,Q]/6p{r) is given by the functional 
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derivative of Ty[p,Q]. It appears because of the different definitions of the exchange- 
correlation energy in the restricted and unrestricted theories. The component of the 
exchange-correlation potential acting on the spin density is given by 

.;r|p.Q|(.) ^ ' J^f}^f ^ 1 + §1^) . (104) 

For the ground-state electron and spin densities, the first term vanishes according to 



Eq. (74), and the above expression for v^^'^[p,Q]{r) reduces to the Euler-Lagrange equa- 



tion for the spin density [cf. Eqs. (90) and (p3l|]. 



In summary, the KS potential in the spin-unrestricted case differs from the spin-restricted 
KS potential v^dplif) by (a) an additional component fx?"[p, Q]{r) acting on the spin 
density, and (b) a correction to the spin-independent potential Vu[p,Q]{r). Thus, starting 
from a spin-restricted reference system with the total electron density p{r) and with a 
spin density that differs from the spin density Qs{r) of the interacting system, the KS 
potential is modified such that its spin density becomes equal to the one of the fully 
interacting system Q{r). To achieve this, the spin potential v^^dp,Q]{r) az has to be 
introduced. However, with the total potential kept fixed, this would lead to a change of 
the total electron density, and in order to keep p(r) unchanged, the correction Vu[p, Q]{f) 
is needed. 



4.4 Spin States in KS-DFT 

So far, we have only considered a spin-state independent theory, i.e., with the exact 
exchange-correlation functionals, the spin-restricted and spin-unrestricted KS-DFT for- 
malism discussed above will lead to the ground-state, irrespective of its spin symmetry. 



As discussed in Section 3.4, targeting the lowest state of a given spin symmetry (i.e., 
with a specific eigenvalue S{S + 1) of S^) requires a spin-state specific Hohenberg-Kohn 



functional -Fhk[p] defined in Eq. (56). 
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In the spin-restricted case, this is formally possible by employing the spin-state indepen- 
dent definition of the noninter acting kinetic energy, which results in a spin-state specific 
exchange-correlation functional. In practice, a different strategy is followed: The nonin- 
teracting reference system is defined such that it is described by a single Slater determi- 
nant with Ms = S. This can be achieved by defining the spin-state specific noninteracting 
kinetic energy as 

Tf [p] = min (^^'^s=s^f^^Ms=s^^ ^^q^-^ 

where $^s=s jg Slater determinant with Ms = S. Such a Slater determinant is always 
an eigenfunction of S'^ with eigenvalue S'(S'+ 1). Thus, it is ensured that the wavefunction 
of the noninteracting reference system has the same spin symmetry as the wavefunction 
of the true interacting system. However, as always in spin-restricted KS-DFT, for S* > 
the spin density Qf^^^^ of the noninteracting reference system will differ from the one of 
the fully interacting system. 

With this definition of a spin-state specific noninteracting kinetic energy, the spin-state 
specific exchange-correlation energy is given by 

E^[p] = Fi^[p]-T^[p]-J[p]. (106) 

One protocol for constructing this spin-state dependent exchange-correlation functional 
then proceeds by using the spin density Qg^"^^^ of the noninteracting reference system to 
distinguish the different spin states. To this end, E^^[p] is expressed as 

E!Ap] = E^c'^[p,Qf'='[p]]- (107) 

Here, Eitf\p, Q] is not equal to the spin-resolved exchange-correlation functional E^dp, Q] 



defined in Eq. (68). Its dependence on Q does not describe the spin-density dependence of 



the exchange-correlation energy, but instead introduces the spin-state dependence. This is 



indicated by the superscript "(ss)" . The multiplet-DFT scheme of Daul 117,118 and the 



restricted open-shell KS (ROKS) scheme 119 124 as well as related approaches 125,126 
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proceed along these lines, but usually include additional ideas originating from Hartree- 



Fock theory 127,128 



It appears that if applied in a spin-restricted formalism, all available approximate exchange- 
correlation functionals have to be understood as approximations to E^c^ [p, Q] and not as 
approximations to Exc[p,Q]- This has important consequences for the construction of 
such approximate exchange-correlation functional. In particular, one has to realize that 
Exc^ [p, Q] does not fulfill the fractional spin condition discussed in Section 



3.3 



and thus 



this condition should not be included when constructing approximations to it. 



In addition, in spin-restricted KS-DFT the ground-state spin density is not directly avail- 
able and has to be determined after calculating the ground-state electron density by 



minimizing E^dPiQ] with respect to Q as discussed in Section 4.1.4 In this step, one 
has to use E^c[p,Q] — which includes the correct spin-density dependence — instead of 
Exc^ [p, Q] ■ Thus, the construction of a different class of approximate exchange-correlation 
functionals that include the spin-density dependence by approximating Exc[p,Q] (or its 
spin-state specific analogue E^^[p, Q]) instead of E^c \p, Q] would be required. Of course, 
the fractional spin condition applies to E^dp, Q] and E^^p, Q], and should also be incor- 
porated when constructing such approximations. 

In spin-unrestricted KS-DFT, it is not possible to require that the noninteracting reference 
system is an eigenfunction of S"^ 88,112 . Thus, it is not possible to define a spin-state 



specific analogue of t]"^ [p, Q] . Consequently, the spin-state dependence only enters the 
exchange-correlation functional, which becomes 

Eif''[p,Q] = Fi^[p,Q] - ri")[p,g] - J[p]. (108) 

For constructing approximations to this spin-state specific exchange-correlation func- 
tional, usually a strategy similar to the one in spin-restricted KS-DFT is applied. To this 
end, the description is restricted to the case of Ms = S, i.e., only the maximal eigenvalue 
of Sz is allowed. Then, the spin-state specific exchange-correlation functional can be 
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expressed as 

E^^''[p,Q] = Eir'^[p,Q''s=sy (109) 

Here, different spin states can be distinguished based on the integral of the spin density. 
However, the functional E^^'^'^[p,Q] does not describe the correct spin-density depen- 
dence of E^\p,Q] anymore. Neither does El^^'^[p,Q], which is limited to spin densities 
corresponding to Ms = S. However, the spin densities Q^^^ corresponding to other eigen- 



values of Sz can be obtained from the scaling relation of Eq. (33). Again, it is important 
to realize that the fractional spin condition does not apply to E^'''^^[p,Q]. 

The idea of using the spin density as a means to distinguish different spin states in 
a spin-unrestricted KS-DFT formalism is taken even further in broken-symmetry DFT 



24 , 40 , 44 , where the requirement that the spin density of the noninteracting reference 
system matches the correct spin density of the fully interacting system is sacrificed in 
favor of obtaining accurate energetics for low-spin states. Consequently, it has been 
suggested that in this case, the spin density in fact serves to describe the (spin-state 



specific) on-top pair density 129 . If broken-symmetry DFT calculations are interpreted 
in this way, one would need to determine the spin density in a separate step from the 
minimization of E^dp, Q] (or its spin-state specific analogue E^^p, Q]), as discussed above 
for spin-restricted KS-DFT. 



5 Spin in Relativistic DFT 

5.1 Spin and Current in Relativistic Quantum Mechanics 

So far, the discussion has focused on spin in nonrelativistic DFT. For the sake of com- 
pleteness, we now consider the generalization to the more fundamental relativistic regime. 



The relativistic theory relies on Dirac's semi-classical theory of the electron (see Ref. 57 
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for a detailed account). This quantum theory describes the relativistic motion of the 
electron in a classical external electromagnetic field, represented by the scalar and vector 
potentials 0cxt(^) and Acxt(^)- The equation of motion reads in this case (in Gaussian 
units), 



h^ij{r,t) 



cct ■p + (3meC^ + gg 



^{r,t)=ih^^^{r,t), (110) 



which yields an energy eigenvalue equation for determining the stationary states when 
the right hand side is replaced by Eilj{r,t). The Dirac Hamiltonian consists of the 
kinetic energy operator ca ■ p, the rest energy term PrrieC^, and the interaction operator 
Qe{(pcxt{i") — CK • A(r)), where c is the speed of light and p = —iKS7 is the momentum 
operator. The Dirac matrices are contained in the parameters ck = [ax,(yy,az) and /3, 
which are in the standard representation. 



a. 



a 



Or 



and 



:iiii 



where cxj are the the Pauli spin matrices ax, Uy, and o"^ defined in Eq. Q. A consequence 
of these four-dimensional operators is that the Dirac Hamiltonian /i^ is a 4 x 4 matrix 
operator with a four-component eigenvector ilj{r,t), the so-called 4-spinor. Such a four- 
component description naturally includes spin, which had to be introduced in an ad hoc 
fashion in the nonrelativistic theory. 



In analogy to the nonrelativistic case, one can define the relativistic spin operator as 



h 



.(4) 



h 



r(4) ^(4) ^(4)aT 



^y^x 



with cr, 



(4) 



CTi 



(112) 



This operators still obey the commutation relations of an angular momentum. However, 
in the relativistic case and do not commute with the Hamiltonian. Therefore, 
eigenstates of cannot be chosen as eigenfunctions of and anymore and spin is 
thus not a good quantum number in relativistic theory. In spherically symmetric systems 
such as atoms, one can consider the total angular momentum instead. This also has the 
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consequence that, in contrast to the nonrelativistic case, the expectation value of will 
depend on the choice of the quantization axis. 

In the case of many electrons, the relativistic wavefunction again assumes a tensor struc- 
ture, i.e., for electrons, the wavefunction formally has 4^ components. It has to fulfill 
the Pauli principle by being antisymmetric with respect to the exchange of any two elec- 
trons. The relativistic many-electron Hamiltonian is then given by (neglecting projectors 
on positive-energy solutions for the sake of brevity), 

N N N 

H'' = Y.^h^^) + Y,Y. (113) 

i=l i=l j=i+l 



where h^{i) is the one-electron Hamiltonian of Eq. (110) acting on electron i and g{i,j) 
is the operator describing the interaction between electrons i and j. The form of the ex- 
act electron-electron interaction operator can be derived from quantum electrodynamics, 
but usually only approximate forms are employed in practice. The simplest approxima- 
tion is to employ the nonrelativistic Coulomb operator. However, the resulting Dirac- 
Coulomb Hamiltonian is not Lorentz invariant. This approximation is improved by the 
Dirac-Coulomb-Gaunt Hamiltonian, which also includes the unretarded magnetic inter- 
action between electrons, whereas the Dirac-Coulomb-Breit Hamiltonian approximately 
describes the retarded electromagnetic interaction (for a detailed discussion, see chapter 8 
in Ref. [57]). 



As in nonrelativistic theory, an electron density and a current density can be defined in 
the relativistic many-electron theory such that they fulfill a continuity equation. This 



results in the definition of the electron density as 57, 130 

p(r) =N j ^^{r,r2,..., r^) ^(r, rs, . . . , r^) dVs ■ ■ ■ d^r^ (114) 
and of the current density as 

j{r) = cN y"^"f(r,r2,...,r^)cKi^(r,r2,...,r^)dV2---dV^ (115) 
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where ai indicates that the Dirac matrices act on the first electron and the dagger denotes 
the transposed and complex conjugate spinor. This definition of the current density still 
holds in the presence of external magnetic fields. 

5.2 Relativistic Current-DFT 

A central aspect of relativistic theories is that all fundamental physical equations must 
preserve their form under Lorcntz transformations from one incrtial frame of reference to 
another one. This requires that the equations are castable in tensorial form. Hence, all 
quantities are joined to 4- vectors, which are basic physical quantities in any relativistic 
theory. For instance, the electromagnetic potentials arc joined to yield the 4-potential 
A^^ — {4>, A). This is the reason why a relativistic theory has to include both the scalar 
and the vector potentials simultaneously. The density p{r) is also part of a 4-vector, 
namely of the 4-current j^{r). The other three components are given by the current 
density j{r) such that — {cp,j). 

Now, a similar decomposition of the energy expectation value as in the nonrelativistic 
case can be performed, 

E = (*|^^|*) = (*|r^|*) + (^iKxtl*) + (*|v;e|*), (116) 

where is the relativistic "kinetic energy" operator, collecting all terms of the one- 
electron Hamiltonian containing the Dirac matrices a and (3 (i.e., the kinetic energy 
and the rest energy terms), Kxt consists of the remaining one-electron terms and describes 
the interaction with the external electromagnetic potentials, and Vge is the electron- 
electron interaction operator. The second term can be calculated directly from the 4- 
current, without the need to know the full wavef unction, as 

(*|V;xt|*) = | j Ur)AUr)d'r 

^q,J p(r)0ext(r-) + | / Ji^) ■ ^ext(r-) dV = Kxtb'1, (117) 
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where Einstein's convention of implicit summation over repeated lower and upper indices 
has been employed. The first term is the same as in the nonrelativistic case and describes 
the electrostatic interaction of the electron density with the external potential, whereas 
the second term accounts for the interaction of the current density with the external 
magnetic field. 

A 4-current may be understood as the source of a 4-potential which can be calculated 
from the relativistic generalization of the Poisson equation of electrostatics, 

UA^{r) = ^f{r), (118) 

with the D'Alembertian operator being the Minkowski space generalization of the three- 
dimensional Laplacian, i.e., it is defined as □ = — A. From the solution of this 

equation, one obtains for the 4-potential generated by the electronic 4-current j^, 

The classical interaction energy Jfj'^] of the 4-current j'^ with itself (i.e., the relativistic 
analogue of the classical Coulomb interaction in nonrelativistic theory) is then given by 

2c^ J \r — r'\ 

= ^ / dVdV + 4 / '^P^ dVdV'. (120) 

2 J \r-r'\ 2c2 J \r - r'\ ^ ^ 

Here, the first term is the classical electrostatic (Coulomb) interaction J[p], whereas the 

second term accounts for the magnetic interaction of the electrons. Note that this second 

term should only be present if the Gaunt or Breit interaction has been included in the 

Hamiltonian. 



The analogy between Kxtb^] and J[j^] [Eqs. (117) and (120)] and their nonrelativistic 
counterparts Kxt[p] and J[p] suggests that the fundamental quantity for a relativistic 
formulation of DFT is actually the 4-current, i.e., the combination of the electronic density 
p(r) and the current density j{r). For this reason, the relativistic generalization of 



standard DFT has been known as current-density functional theory (CDFT) [131 134 
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It can be shown that in relativistic theory, an analogue of the HK theorem exists: The 
4-current j'^(r') uniquely determines the external 4-potential A^{r) up to a gauge trans- 
formation |131| . Therefore, it is possible to consider the energy E[j^] as a functional of 
the 4-current j^{r). As in the nonrelativistic case, one then proceeds by introducing a 
reference system of noninteracting electrons with the relativistic Hamiltonian 

N 

Hf = Y,cai^)■Pi + Pi^^eC^ + (le[Mr^)-a{^)■A,{ri)y (121) 

i=l 

The eigenfunctions of this noninteracting Hamiltonian are given by Slater determinants 
"^s = \ fi,f2, ■ ■ ■ ,V^iv| constructed from 4-spinors '^i{r), which solve the relativistic one- 



electron KS equations 



ca p + jSnieC^ + Qe {4>s{r) - OL ■ As{r] 



^,{r) = eiifiiir). (122) 



For the resulting single Slater determinant, the electron density is given by 

occ 

p{r) = J2^K^)v^{r), (123) 

i 

an the relativistic current density reads 

occ 

j{r) = cJ2'PKf)^'P^if)■ (124) 

i 

The noninteracting kinetic energy Ts[j'^] can then be defined as the kinetic energy = 
(^s|T-^|\E's) of such a noninteracting system with the 4-current j'^, which allows for a 
similar decomposition of the total energy functional as in nonrelativistic DFT, i.e., 

E[f] = T,[r] + Kxt[j^] + m + i^xcb'l- (125) 

Here, the exchange-correlation energy functional is defined as containing all the energy 
contributions not accounted for by the first three terms. In analogy to nonrelativistic 
KS-DFT, it can then be shown that the 4-current j^(r) of the true interacting system 
can be determined from the self-consistent solution of the relativistic KS equations for a 



noninteracting system 131 135 , with the 4-potential 



AW = AUr) + A^[j^]ir) + (126) 
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The functional derivatives of the exchange-correlation functional with respect to the com- 
ponents of the 4-current defines the exchange-correlation potential v^^[j^]{r), which now 
assumes a four-component form. 



5.3 Electron Density and Spin Density in Relativistic DFT 

In view of the previous discussions concerning nonrelativistic KS-DFT, we now face the 
following question: (i) Is it still possible to formulate a relativistic DFT in terms of the 
electron density only, (ii) how is the fundamental 4-current related to the density and 
spin density considered as fundamental quantities in nonrelativistic DFT, and (iii) how 
do the nonrelativistic restricted and unrestricted formulations of KS-DFT emerge from 
the relativistic framework? 

Regarding the first question, under some additional assumptions it is indeed possible 



to prove a relativistic HK theorem for the electron density only 135 137 : In the case 
considered throughout this paper, the external potential is the electrostatic potential of 
all atomic nuclei in a molecule that are at rest (Born-Oppenheimer approximation). As 
the nuclei are not moving, they do not create magnetic fields. In addition, we neglect 
any magnetic fields that stem from nuclear spins and assume that there are no addi- 
tional external electromagnetic fields. As a consequence, the external electromagnetic 
4-potential only contains the time-independent scalar potential of the atomic nuclei, 
i.e., (pextif) = Vnuc{r) and A{r) = 0. This assumption is also common practice in almost 
every relativistic quantum chemical calculation. Then, it can be shown that the external 
scalar potential 0ext(^) within this specific reference frame is (up to a constant) uniquely 



determined by the electron density p(r) only 135 137 . Within such a framework, the 
relativistic total energy functional becomes 

E[p] = T^lflp]] + KxtM + Jinp]] + E^^irip]]. (127) 
Here, the current density j{r) is still required, but it is now uniquely determined by the 
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electron density (i.e., j = j[p]), in the same way in which the spin density is determined 
by the electron density in nonrelativistic DFT. If a magnetic interaction between the 
electron is not included (i.e., the Dirac-Coulomb Hamiltonian is employed), also the 
analogue of the Coulomb interaction reduces to a functional J[p] of the density only |135|. 
Moreover, it now becomes possible to set up theories in which not the full 4-current is 
used as fundamental variable, but only the parts of it that are related to the total electron 
density and the spin density. This will be discussed further in the following subsection. 

In order to understand the relation of the 4-current to the spin density, it is important to 
realize that the definition of the current density (naturally) involves a velocity operator, 
which is in close analogy to classical mechanics (correspondence principle) [57]. The 
velocity operator in Dirac's theory of the electron follows from the Heisenberg equation 
of motion applied to the position operator and turns out to be ccx.. Hence, as the Dirac a 
matrices contain the Pauli spin matrices cr, we see immediately that the current density 



j carries the spin information 130 



The relation to the spin density can be made more explicit by invoking a Gordon com- 
position of the current density which separates it into a charge- and a spin-related cur- 
rent 



53 ,131 . For the one-electron case, one can carry out this decomposition by rewriting 



the Dirac eigenvalue equation as 
1 



^(r) 



:i28) 



and by splitting up the definition of the current in a somewhat artificial way as 



c c 

- ip^r) Cfip{r) + - ip^r) cx'ip{r). 



(129) 



Eq. (128) can now be used to replace ip in the first term and ip"^ in the last term. As is 



shown in the Appendix, by exploiting the commutation relations of the Dirac matrices 
{a'', /9} = and a^a'- = Ma^, a'] -|- Ma^, a'} = iekinO'n^ + ^kh one arrives at 



53 



131 



see 
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also pages 552-558 in Ref. 138| ) 

. 1 



2mp 



(130) 



where the 3- vector cr^^^ = (ai , cr- 



(4) ^(4) 



contains the 4x4 Pauh matrices a,-^'* introduced 



in Eq. (112). In a many-electron system described by the Dirac-Coulomb Hamiltonian, a 



similar decomposition of the current density can be performed (even though the derivation 
becomes slightly more complicated, see Appendix), and one obtains. 



J(nj 



J ^t/3icrf)^dV2---dV. 



dVa ■ ■ ■ d^Tv 
(131) 



The first line of this expression resembles the definition of the current density in nonrela- 
tivistic quantum mechanics, whereas the second term can be identified as arising from the 
electron spin. This can be made more apparent by defining the magnetization (density). 



m(ri) = N j¥ ^ ^ dVa ■ ■ ■ d^ 



(132) 



Then, the contribution of the second term of Eq. (131 ) to the interaction energy with the 



external electromagnetic potentials [cf. Eq. (117)] becomes 

V:r[f] = ^ / (V X m(r)) ■ Ae.t(r) dV = -//^ J m(r) ■ B,^,{r) d\ (133) 

where Bext = V x Acxt is the external magnetic field. The minus sign originates from 
the negative charge of the electron. This closely resembles the form of the spin Zeeman 



interaction in the nonrelativistic case [cf. Eq. (30)]. If an inhomogeneous magnetic field 



in z-direction is considered, Eq. (133) reduces to 



(134) 



and by comparison with Eq. (31) we notice that rrizir) can be identified with the spin 



density Q{r) in nonrelativistic theory. However, while in the nonrelativistic case the spin 



50 



density is (in the absence of external magnetic fields) independent of the choice of the 
quantization axis, this is not the case in the relativistic theory, where does not commute 
with the Hamiltonian because of the presence of spin-orbit interactions. 

5.4 Relativistic Spin-DFT 

In the relativistic CDFT formalism discussed above, the noninteracting reference system 
is chosen such that it has the same 4-current j^(r') [i.e., the same electron density p(r) 
and the same current density jir^ as the true interacting system. This results in a 
noninteracting kinetic-energy functional T's[p,j] that can be defined as 

T,[p,i]= min (^.If^l^',), (135) 

where the constrained search has to be restricted to positive-energy wavefunctions ^5 
to avoid a variational coUaps. In such a formalism, the KS equations then contain a 
four-component exchange-correlation potential. 

However, we also pointed out that in the case of molecular systems in the absence of 
external magnetic fields, a description relying on the electron density p{r) only as funda- 
mental variable is sufficient. Therefore, it is formally also possible to develop a relativistic 
"density-only" KS-DFT that resembles the nonrelativistic restricted KS-DFT formalism. 
This can be achieved by only requiring from the noninteracting reference system that it 
has the same electron density as the interacting system, and consequently defining the 
noninteracting kinetic energy as, 

Tf\p\ = min (*i'^)|r|*('^)>, (136) 

again restricting ^ g to positive energy solutions. Then, the total energy functional can 
be decomposed as, 

m = Tf^ \p\ + Kxt \p\ + J\P\ + [P] , (137) 
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where we have neglected the magnetic interactions in J[p\. Note that such a decom- 
position imphes a different definition of the exchange-correlation energy. The result- 
ing KS equations then feature a one-component exchange-correlation potential v^c [p\ = 
{1/ qe) 6E^[p\/ 6p{r). However, in such a formalism the current density 3s{r) and the 
magnetization ms{r) of the noninteracting reference system do not agree with the true 
interacting system. Instead, these are again a functional of the density only, just as the 
spin density is in nonrelativistic restricted KS-DFT. 

In between full CDFT and relativistic density-only KS-DFT, different intermediate for- 
mulations of relativistic KS-DFT are now also possible. For a relativistic system of non- 
interacting electrons, the magnetization is given by 

occ 

m{r) = Y,v\{r)PcT^'^^^{r). (138) 

i 

The noninteracting reference system can then be set up such that, in addition to the elec- 
tron density, some parts of the magnetization density agree with those of the interacting 
system. For instance, we can require that the z-components of the magnetization m^(r) 
match, or we can demand that the lengths of the magnetization vector |m(r)| at each 
point in space agree. The choice of what quantity is to be reproduced by a relativistic 
KS-DFT formalism in addition to the density is our freedom of choice. If we choose to 
also reproduce the 2;-component of the magnetization density, this approach is called the 
collinear ["(c/)"] approach because an artificial external global quantization axis is intro- 
duced for the spin. If we require to reproduce the length of magnetization instead, this is 
a noncollinear ["(nc)"] approach because the magnetization is a vector field representing 



a magnetic dipole moment whose direction changes with position 135, 137, 139 



141 



Note that these different choices for the noninteracting reference system each implies 
a different definition of the noninteracting kinetic energy functional, Ti^''' [p, m^] and 
Ti"^-* [p, |m|] , and thus also of the exchange-correlation energy, E^P [p, and -EI"^'' [p, \m\\, 
respectively. For the choices mentioned here, where in addition to the electron density one 
additional quantity is reproduced by the KS system, the resulting exchange-correlation 
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potential has two components. This is in close analogy to the case of nonrelativistic unre- 
stricted KS-DFT. Therefore, approximate exchange-correlation functionals developed in 
the nonrelativistic domain are usually employed in practical applications of such relativis- 
tic spin-DFT schemes. However, the exchange-correlation potential is defined differently 
in nonrelativistic unrestricted KS-DFT and in the relativistic collinear and noncollinear 
cases. Consequently, also different exact conditions apply to the exchange-correlation 
functional. The different possible choices for setting up relativistic KS-DFT are summa- 
rized in Fig. |3] 

Finally, we discuss how the nonrelativistic unrestricted KS-DFT formalism emerges from 
relativistic spin-DFT. This is most easily seen for the collinear approach, although the 
noncollinear one reduces to the same nonrelativistic limit as well. The 2;-component of 
the magnetization density of the KS reference system is given by 



{r) = J2vl{r)fiai'^^.{r). (139) 



If spin-orbit coupling is neglected, the Hamiltonian commutes with the spin operators 
and the KS spinors can each be expressed as (two-component) spin orbitals ipi^a = f^ci 
or (fi^j3 = (fiP, where y?"" are 2-spinors consisting of an upper and a lower component. 
Then, the spin-orbit coupling free (SOfree) z-component of the magnetization resembles 
the nonrelativistic spin density, 

occ a P 

- E - \^r''ir)\'] - E [br(r)P - |^f"(r)p] = Q(r), 

i i 

(140) 

where the superscript "f/" and "L" denote the upper and lower components, respectively. 
The neglect of spin-orbit coupling is, of course, an approximation that yields a scalar- 
relativistic Hamiltonian (which considers kinematic relativistic effects only) or even a 
nonrelativistic Hamiltonian if the speed of light is taken to be infinity. 

53 



The two 2-spinors and Lp^ considered here can each be reduced to one-component spin- 



orbitals if a unitary transformation 142 is performed to decouple the upper and lower 
components. This then also requires a unitary transformation of the operators involved 



in the calculation of 143 . Still, taking the limit c — ?■ oo yields the nonrelativistic 
theory and relativistic spin-DFT reduces to nonrelativistic unrestricted KS-DFT. 



6 Future Directions for Spin-DFT 

As we have seen, different options exist for setting up KS-DFT for open-shell systems. In 
the nonrelativistic case, one has to choose between a spin-restricted and a spin-unrestricted 
formulation of KS-DFT. In the former case, the wavefunction of the noninteracting ref- 
erence system can always be chosen as an eigenfunction of S^, but its spin density 
differs from the correct one. Alternatively, in spin-unrestricted KS-DFT the noninter- 
acting reference system has the correct spin density, but cannot be an eigenfunction of 
S"^. In relativistic DFT even more options in between density-only KS-DFT (with a 
one-component exchange-correlation potential) and full CDFT (with a four-component 
exchange-correlation potential) are possible. In particular, one can choose to reproduce ei- 
ther a single component or the magnitude of the magnetization vector in collinear and non- 
coUinear relativistic KS-DFT (which both employ a two-component exchange-correlation 
potential), respectively. 

Which choices we make may be determined by our ability to set up a proper approx- 
imation to the corresponding exchange-correlation energy functional. However, it is 
important to understand that these different choices imply different definitions of the 
noninteracting kinetic energy and the exchange-correlation functional. Thus, different 
exact conditions apply to these functional, which must be considered when developing 
such approximations. This is most obvious for the fractional spin condition in the nonrel- 
ativistic case. While in spin-restricted KS-DFT, the constancy condition directly applies 
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to the exchange-correlation functional, in the spin-unrestricted case it does not hold for 
the exchange-correlation functional alone, but to the sum of noninteracting kinetic and 
exchange-correlation energy. It can be expected that in the latter case, devising approx- 
imate exchange-correlation functionals that include this condition will be significantly 
more difficult. 

While exact DFT should always lead to the correct ground-state — irrespective of its 
spin state — in practice it appears more useful to rely on a theory that is able to target 
different spin states separately. This should also simplify the development of approximate 
functionals because it becomes possible to account for the different exact conditions ap- 



plying to the exchange-correlation hole for different spin states 129 , 144 . Note that such 
a spin-state specific DFT always has to operate within a spin-orbit coupling free frame- 
work, because in the fully relativistic theory S"^ does not commute with the Hamiltonian 
and different spin states do not correspond to the lowest state of a specific symmetry 
anymore. 

Thus, besides finding approximations that accurately account for the spin-density de- 
pendence of the (nonrelativistic) exchange-correlation functionals E^c[p,Q] or E^\p,Q], 
finding ways of including the spin-state dependence into these functionals is another im- 
portant problem for open-shell systems. In common approximations, it appears that the 
spin-density dependence is actually used to model this spin-state dependence (i.e., the 
integral of the spin density is used to distinguish spin states). Spin-state and spin-density 
dependence of the exchange-correlation functional are intermingled in all available ap- 
proximate functionals, which manifests itself in their violation of the fractional spin con- 



dition 101 , 102 . To make progress in the development of reliable density-functional ap- 
proximations, we believe it will be essential to consider the spin-state and the spin-density 
dependence of the exchange-correlation functional separately. 
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A Gordon Decomposition of the Current Density 



In the one-electron case, the current density is given by 



j(r) = cilj{ryailj{r). (141) 
For the fc-component, we can rewrite this definition in a somewhat artificial way as 

/ = cilj^a'^ilj = - V^^aV + - i^^a'^i/j, (142) 

where we dropped the dependence of on the spatial coordinate to simplify the notation. 
The eigenvalue equation of the one-electron Dirac Hamiltonian equation can be rewritten 
to obtain an expression for ip, 

1 



(143) 



and by taking the transpose and complex conjugate also an expression for ip'^ (exploiting 
1 



t 



(144) 



These expressions can now be used to replace ip in the first term and ip'^ in the second 



term of Eq. (|142|) to obtain, 

c(3cx- {p-'^-^A) + P{e -q.cp) 
^p* - '^A}j^^ ■ot(3+[E- qe(p)^^(3 



2meC 
1 



2mpC 



(145) 
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and by reordering the different terms we get 



1 



2mp 



+ 



2mpC ( ^^^) 



(146) 



The last term is zero because of {a'^, ^} = 0. For the first term, we use 



(147) 



to arrive at (employing the convention of implicit summation over repeated indices) 
1 



2me 
1 

2mp 



After regrouping the different terms, we obtain 
1 



2mp 



+ 



+ 



1 



2mp 



2mp 



(4)/_* 



where the last term is zero because Skin = —^ikn- Hence, we find 



2mp 



+ 



2mp 



+ 



d 



2me 9n 

and combining the different components of j again gives 



2mp 



h 



(148) 



(149) 



(150) 



+ -^V X V'V/3(t(^))V^. (151) 
zrup 
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For an A^-electron system, we start from the definition of the many-electron current density 



130 



j(ri) = j ^(ri, rs, . . . , cki ^(ri, rs, . . . , vn) dVs ■ ■ ■ d' 



(152) 



and the wavefunction is an eigenfunction of the many-electron Dirac-Coulomb Hamilto- 



nian 



TV 



(153) 



1=2 



projected onto the electronic (positive-energy) states. 



We can now rewrite the corresponding eigenvalue equation in a similar fashion as in the 
one-electron case to obtain for the wavefunction 

i=2 



N 



-choLx- (pi - |^(n)) (i?-ge0(ri)) /i^(z) -/^iK 

and for its transpose and complex conjugate 

i=2 



m (154) 



vl/t 



pt-|A(ri))^t«i/3i+(i?-g,0(ri))^t/3i-5^ (/i^(z)vl/) Vi-(Ke^)^/3 

(155) 



Now, we can rewrite the fc-component of the integrand in Eq. (152) as 



(156) 



After substituting Eq. (154) for in the first term and Eq. (155) for in the second 



term, we obtain, in addition to the term already present in the one-electron case [cf. 



Eq. 1^], 



1 



N 



2mpC 



vi/t«tLgi^(/i^(^)vi>)^ + /3iV;e^ 



i=2 



1 ^ 

i=2 



N 



N 



«=2 



i=2 



2meC . 
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With a multiplicative operator acting equally on all components of the wavefunction for 
Vee, the last term is zero because of {a^, 13} = 0. However, if the Gaunt or Breit interaction 
is included, the corresponding term is not zero and gives an additional contribution to 
the current that is not present in the one-electron case. This term contains the spin-spin 
interactions. 

For the remaining terms, we now consider one of the terms separately, for instance h^{'2), 
and note that 



2meC 



2nipC 



(157) 



because of {aj^,/3i} — 0. Integrating the expression in square brackets over r2 we now 
find 

^^at/3i(/i^(2)^) d^r2- J (/i^(2)^)^/3i^ d^s 
*ta^/3i(/i^(2)*) dVs- / ^^aj^/3i(/i^(2)^) dVs = 



... - (158) 

where for the second term we exploited that h^(2) is hermitian and commutes with 
and Pi because these act on different electrons. Combining all of these results, we find 
that (with the Dirac-Coulomb Hamiltonian) , 

j{ri)^cN J *(ri,r2, . . . ,rjv)^ai*(ri,r2, . . . ,rjv) dV2---dVjv 

^^(3i{pi - |A(ri))* + (pl - |A(ri))*t^i* 



2mp 



N 



+ AT / V X (^t p^f^W ^1^2 • • • dVjv 
2me J 



d^r2 • • • d^rjv 
(159) 
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Figure 1: Connection between the electron density p{r) and the total energy as stated by 
the first Hohenberg-Kohn theorem. 

a) Hohenberg-Kohn theorem for i;-representable densities 
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b) Hohenberg-Kohn theorem for A^-representable densities 
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Figure 2: Relationship between wave function, total electron density, and spin density of 
the system of fully interacting electrons and of the spin-restricted and spin-unrestricted 
Kohn-Sham reference systems of noninter acting electrons. 
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Figure 3: Relationship between wave function, total electron density, and the magnetiza- 
tion of the system of fully interacting electrons and of the Kohn-Sham reference systems 
of noninteracting electrons in different version of relativistic KS-DFT. 
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Table I: Comparison of the spin-restricted and spin-unrestricted formulations of KS-DFT. 
"Correct" indicates that the quantity calculated for the noninteracting reference system 
agrees with the corresponding one of the fully interacting system. 
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